Each of the models in shadia can be run under a “no dams” scenario to understand the operational baselines for each population included in the package. In other words, by specifying upstream passage probabilities of 1.00 at all dams in a given river, the user can simulate a scenario in which there is unimpeded passage to spawning habitat.
This is easily achieved for any of the models by running with the default parameter values.
Here is an example using the Penobscot River, with all other user-specified arguments set to their default values:
We could then get the output and plot our result. Of course, this is only a single run, so your results may be very different from what is shown here.
result = nodams$res
plot(x = result$year, y = result$populationSize,
type = 'l', lwd=2, col = 'red', xlab='Year',
ylab='Spawner abundance')The above example would run the model a single time for the default number of years (50, in the case of penobscotRiverModel()). If we wanted to do this several hundred times to get outputs that are robust to variability in the input distributions, then we could wrap the call above in a snowfall script.
This might look something like:
# Load R packages
library(snowfall)
library(rlecuyer)
library(shadia)
library(plyr)
# Initialize parallel mode using sockets
sfInit(parallel=TRUE, cpus=7, type="SOCK")
wrapper <- function(idx) {
# Run the model
res1 <- penobscotRiverModel(
nRuns = 1,
nYears = 50,
timing = list(1,1,1,1,1,1,1),
upstream = list(
milford = 1,
howland = 1,
westEnfield = 1,
brownsMill = 1,
moosehead = 1,
guilford = 1,
weldon = 1
),
downstream = list(
stillwater = 1,
orono = 1,
milford = 1,
howland = 1,
westEnfield = 1,
brownsMill = 1,
moosehead = 1,
guilford = 1,
weldon = 1
),
pinHarvest = 0,
inRiverF = 0,
commercialF = 0,
bycatchF = 0,
indirect = 1,
latent = 1,
watershed = TRUE
)
# Define the output lists
retlist <- list(sim=res1)
return(retlist)
}
# Export needed data to workers
# load required packages on workers.
sfLibrary(shadia)
# Distribute calculation to workers
niterations <- 30
start <- Sys.time()
# Use sfLapply() function to send wrapper() to the workers:
result <- sfLapply(1:niterations, wrapper)
# Stop snowfall
Sys.time()-start
sfStop()
# Extract results list from output list
out <- lapply(result, function(x) x[[c('sim')]])
# Extract user inputs and population metrics
res <- lapply(out, function(x) x[[c('res')]])
resdf <- do.call(rbind, res)
# Extract sensitivity variables
sens <- lapply(out, function(x) x[[c('sens')]])
sensdf <- do.call(rbind, sens)
# Have a look at result
plotter <- ddply(resdf, 'year', summarize,
mu=mean(populationSize),
lwr=quantile(populationSize, probs=0.25),
upr=quantile(populationSize, probs=0.75)
)
plot(plotter$year, plotter$mu, type = 'l',
xlab = 'Year', ylab='Spawner abundance',
ylim=c(0, max(plotter$upr)))
polygon(x=c(plotter$year, rev(plotter$year)),
y=c(plotter$lwr, rev(plotter$upr)),
col='gray87', border=NA)
lines(plotter$year, plotter$mu, col='red', lwd=2)This work is licensed under a GNU General Public License v 3.0.